home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / cpp_libs / rwvector.lha / RWVector2.1 / src / xludecmp.cc < prev    next >
C/C++ Source or Header  |  1989-08-18  |  3KB  |  121 lines
  1. /*
  2.  *    <T> Precision LU Decomposition
  3.  *
  4.  *    Copyright (C) 1988, 1989.
  5.  *
  6.  *    Dr. Thomas Keffer
  7.  *    Rogue Wave Associates
  8.  *    P.O. Box 85341
  9.  *    Seattle WA 98145-1341
  10.  *
  11.  *    Permission to use, copy, modify, and distribute this
  12.  *    software and its documentation for any purpose and
  13.  *    without fee is hereby granted, provided that the
  14.  *    above copyright notice appear in all copies and that
  15.  *    both that copyright notice and this permission notice
  16.  *    appear in supporting documentation.
  17.  *    
  18.  *    This software is provided "as is" without any
  19.  *    expressed or implied warranty.
  20.  *
  21.  *
  22.  *    @(#)xludecmp.cc    2.1    8/18/89
  23.  */
  24.  
  25. /*
  26.  *    These routines call the appropriate LINPACK routines
  27.  */
  28.  
  29. #define NO_VECTOR_MATHFUN
  30. #include "rw/<A>LUDecomp.h"
  31. #include "rw/linpack.h"
  32. #include <stdio.h>
  33.  
  34. // Construct an LU decomposition from m
  35. <A>LUDecomp::<A>LUDecomp(const <A>GEMatrix& m)
  36.      : (m.copy())
  37. {
  38.   assertSquare();
  39.   order = m.rows();
  40.   ipvt = new fortran_int[order];
  41.  
  42.   <A>gefa_(data(), &order, &order, ipvt, &info);
  43. }
  44.  
  45. // Get the determinant from an existing LU decomposition
  46. <T>
  47. determinant(const <A>LUDecomp& d)
  48. {
  49.   d.assertDefined();
  50.   d.assertPivots();
  51.   
  52.   <T>        determine[2];
  53.   <T>*        temp= new <T>[d.order];
  54.   fortran_int    job = 10;    // Indicate only determinant to be computed
  55.   
  56.   <A>gedi_(d.data(), &d.order, &d.order, d.ipvt, determine,
  57.      temp, &job);
  58.   
  59.   delete temp;
  60.   return determine[1] * pow(10.0, determine[2]);
  61. }
  62.  
  63. // Calculate the inverse of an existing LU decomposition
  64. <A>GEMatrix
  65. inverse(const <A>LUDecomp& d)
  66. {
  67.   d.assertDefined();
  68.   d.assertPivots();
  69.   
  70.   <T>        determine[2];
  71.   <T>*        temp= new <T>[d.order];
  72.   fortran_int    job = 1;    // Indicate only inverse to be computed
  73.   <A>GEMatrix    the_inverse = d.copy();
  74.   
  75.   <A>gedi_(the_inverse.data(), &d.order, &d.order, d.ipvt, determine,
  76.      temp, &job);
  77.   
  78.   delete temp;
  79.   return the_inverse;
  80. }
  81.  
  82. <T>Vec
  83. solve(const <A>LUDecomp& d, const <T>Vec& rhs)
  84. {
  85.   d.assertDefined();
  86.   d.assertPivots();
  87.  
  88.   <T>Vec temp = rhs.copy();
  89.   fortran_int job = 0;
  90.  
  91.   <A>gesl_(d.data(), &d.order, &d.order, d.ipvt, temp.data(), &job);
  92.  
  93.   return temp;
  94. }
  95.  
  96. /************ Utility routines ******************/
  97.  
  98. void
  99. <A>LUDecomp::assertDefined()
  100. {
  101.   // Check to make sure this isn't a null matrix
  102.   if(!order){
  103.     char msg[40];
  104.     sprintf(msg, "Undefined LU Decomposition");
  105.     RWnote("<A>LUDecomp::assertDefined()", msg);
  106.     RWerror(DEFAULT);
  107.   }
  108. }
  109.  
  110. void
  111. <A>LUDecomp::assertPivots()
  112. {
  113.   // Check to make sure none of the pivots are zero
  114.   if(info){
  115.     char msg[120];
  116.     sprintf(msg, "Pivot number %d is zero", info);
  117.     RWnote("<A>LUDecomp::assertPivot()", msg);
  118.     RWerror(DEFAULT);
  119.   }
  120. }
  121.